#ifndef am_QSILISTON_HPP
#define am_QSILISTON_HPP


#include <string>
#include <math.h>
#include <tbb\concurrent_vector.h>
#include <tbb\blocked_range2d.h>
#include <boost\date_time.hpp>

#include <borealis\Raster.hpp>
#include "..\AuroraTypes.hpp"


namespace Aurora
{
	/*
	Class: am_QsiListon
		Calculates the incoming shortwave radiation for a given time step following Liston and elder  2002.
		The equation uses two seperate variables depending on the location of the domain. Inland and Costal are the two choices 
		Latitude calculations are currently HARDCODED to Marmot Creek, Alberta. This WILL BE CHANGED IN THE FUTURE.

	In Variables:
		Borealis::Raster<float> Gridded DEM  (meters)
		Borealis::Raster<float> Gridded Slope (degrees)
		Borealis::Raster<float>	Gridded Aspect (degrees, north = 0/360)
		Borealis::Raster<float> Gridded RH
		string					Date header name


	Out Variables:
		Borealis::Raster<float> Qsi

	References:
		Liston, Glen E., and Kelly Elder. 2006. A meteorological distribution system for high-resolution terrestrial modeling (MicroMet). Journal of hydrometeorology 7: 217-234.
	*/
	class am_QsiListon : public AuroraTypes::ModuleBase
	{
	private:

		//in
		int DEM;
		int SLOPE;
		int ASPECT;
		int RH;
		int DATE;


		//out
		int QSI;

	public:

		am_QsiListon(std::string n)
		{
			DEM			= 0;
			SLOPE		= 1;
			ASPECT		= 2;
			RH			= 3;
			DATE		= 4;

			QSI			= 0;

			name = n;
		}

		class CalcQsi
		{
			

			double const m_cosZ;
			double const m_sinZ;
			double const m_Mu;


			boost::shared_ptr<Borealis::Raster<float> > const m_Rh;
			boost::shared_ptr<Borealis::Raster<float> > const m_Qsi;
			boost::shared_ptr<Borealis::Raster<float> > const m_Slope;
			boost::shared_ptr<Borealis::Raster<float> > const m_Aspect;



		public:
			CalcQsi(double cosZ, double sinZ, double mu, boost::shared_ptr<Borealis::Raster<float> > Qsi, boost::shared_ptr<Borealis::Raster<float> > Slope,  boost::shared_ptr<Borealis::Raster<float> > Aspect, boost::shared_ptr<Borealis::Raster<float> > RH)
				: m_cosZ(cosZ), m_sinZ(sinZ),m_Mu(mu),  m_Qsi(Qsi), m_Slope(Slope), m_Aspect(Aspect), m_Rh(RH)
			{

			}

			void operator()(const tbb::blocked_range2d<unsigned int>& r) const
			{

				double cosZ = m_cosZ;
				double sinZ = m_sinZ;

				double mu = m_Mu;

				boost::shared_ptr<Borealis::Raster<float> >  Qsi = m_Qsi;
				boost::shared_ptr<Borealis::Raster<float> >  Slope = m_Slope;
				boost::shared_ptr<Borealis::Raster<float> >  Aspect = m_Aspect;
				boost::shared_ptr<Borealis::Raster<float> >  Rh		= m_Rh;

				
				for(unsigned int x = r.rows().begin(); x!=r.rows().end();++x)
				{
					for(unsigned int y = r.cols().begin(); y!=r.cols().end(); ++y)
					{
						//transmittance
						double slope = Slope->GetRasterAt(x,y);

						//conversion to radian
						slope *= 0.0174532925;

						double aspect = Aspect->GetRasterAt(x,y);
						
						//although north is 0/360, the equations think  south = 0, so add 360/2=180 to it.
						if (aspect >= 180.0) 
							aspect = aspect - 180.0;
						else
							aspect = aspect + 180.0;
						
						aspect *= 0.0174532925;

						double cosi = ( cos(slope) * cosZ + sin(slope) * sinZ * cos(mu-aspect) );
						if(cosi < 0.0)
							cosi = 0.0;
						if(cosZ <= 0.0)
							cosi = 0.0;

						//W/m^2
						//solar extraterrestrial potential
						double Sprime = 1370.0;

						// do the RH conversion


						//cloud cover fraction
						//NEEDS FIXING!!!
						double oc = 0.832*exp((Rh->GetRasterAt(x,y) - 100.0)/41.6);

						double dir = (0.6-0.2*cosZ)*(1.0-oc);
						double dif = (0.3-0.1*cosZ)*oc;

						double qsi = Sprime * (dir * cosi + dif*cosZ);

						Qsi->SetRasterAt(x,y,(float)qsi); 
					}
				}
			}
		};

		void run(AuroraTypes::StationList* stations, AuroraTypes::VariableList* variables, AuroraTypes::InVariables& inVars, AuroraTypes::OutVariables& outVars)
		{

			boost::shared_ptr<Borealis::Raster<float> > Qsi;
			boost::shared_ptr<Borealis::Raster<float> > Dem;
			boost::shared_ptr<Borealis::Raster<float> > Slope;
			boost::shared_ptr<Borealis::Raster<float> > Aspect;
			boost::shared_ptr<Borealis::Raster<float> > Rh;


			try
			{
				Dem = GetRaster(variables,inVars[DEM]);
				Slope  = GetRaster(variables,inVars[SLOPE]);
				Aspect = GetRaster(variables, inVars[ASPECT]);
				Rh = GetRaster(variables,inVars[RH]);


				Qsi.reset(new Borealis::AsciiRaster(Dem->GetRow(),Dem->GetCol()));
				Qsi->SetXllCorner(Dem->GetXllCorner());
				Qsi->SetYllCorner(Dem->GetYllCorner());
				Qsi->SetCellSize(Dem->GetCellSize());


			}
			catch (std::runtime_error e)
			{
				throw std::runtime_error("In " + std::string(__FILE__) +" @ line " + boost::lexical_cast<std::string>(__LINE__) + e.what());
			}


			//pi constant
			double PI = 3.14159;

			//latitude in radians of the domain
			//50 N
			double phi = 0.889364355;

			//day number
			boost::gregorian::date date = stations->at(0)->Now().GetGregorian();

			double d = date.day_of_year();

			/*double hour = to_tm(stations->at(0)->Now().GetPosix()).tm_hour;*/

			double min = to_tm(stations->at(0)->Now().GetPosix()).tm_min;
			double hour =to_tm(stations->at(0)->Now().GetPosix()).tm_hour;

			//latitude of the tropic of cancer
			double phiT = 0.401721464;

			//these are the 2009 values for the summer solstice
			boost::gregorian::date summerSolstice(2009,boost::date_time::Jun,21);
			//summer solstice day of year
			double dr = summerSolstice.day_of_year();

			//solar declination in radians
			double delta = phiT * cos(2.0*PI*( (d - dr)/365.25));

			//hour angle measured from local solar noon
			//double tau = PI *( (hour / 12) - 1) *PI/180.0;
			double tau = (hour * 15.0 -180.0 )*(PI/180.0);

			//Solar zenith angle from local solar noon
			//double Z = acos( sin(delta)*sin(phi) + cos(delta)*cos(phi)*cos(tau) );
			double cosZ = ( sin(delta)*sin(phi) + cos(delta)*cos(phi)*cos(tau) );
			cosZ = max(0.0,cosZ);

			double sinZ = sqrt(1.0 - cosZ*cosZ);

			//solar azimuth
			double mu = asin(max(-1.0,min(1.0,cos(delta)*sin(tau)/sinZ)));
			//double mu = asin( cos(delta)* sin(tau) / sin(Z) );

			//Make the corrections so that the angles below the local horizon
			//  are still measured from the normal to the slope.
			if (tau < 0.0) 
			{
				if (tau < mu) 
				{
						mu = - PI - mu;
				}
				else if (tau > 0.0) 
				{
					if (tau > mu) 
						mu = PI - mu;
				}
			}
						
	

			tbb::blocked_range2d<unsigned int> range(0,Dem->GetRow(),0,Dem->GetCol());
			tbb::parallel_for(range,CalcQsi(cosZ,sinZ,mu,Qsi,Slope,Aspect, Rh));

			AuroraTypes::VariableList::accessor a;
			variables->insert(a,outVars[QSI]);		
			a->second = Qsi;
		}



	};

}



#endif